## This file creates all figures in the main text and appendix ##
## Created by Meredith Dost and last run 8/20/2025 ##

#### SET-UP ####
# load in packages
library(ggplot2)
library(ggpattern)
library(usmap)
library(ggthemes)
library(stringr)
library(dplyr)
library(broom)
library(dotwhisker)
library(cowplot)
library(gridExtra)

# create formatting function for use in plots
dropLeadingZero <- function(l){
  str_replace(l, '0(?=.)', '')
  }

# set working directory
# setwd("")

# load in data--all U.S. counties
burd <- read.csv("final_dataset_for_analysis.csv")
## all border counties w/in 100 miles of border ##
burd100 <- subset(burd, MINDIST<=100)
# read in Medicaid burden measures
kmerg <- read.csv("burden_data/medburden_measures.csv")
# merge in expansion status to Medicaid burden dataset
exp <- read.csv("other_data/Medicaid_expansion_status.csv")
#2014
exp14 <- exp
exp14$exp <- ifelse(exp14$year_implemented==2014, 1, 0)
exp14$exp[is.na(exp14$year_implemented)] <- 0
exp14$year <- 2014
exp14 <- exp14[c("State","year","exp")]
#2016
exp16 <- exp
exp16$exp <- ifelse(exp16$year_implemented>=2014 & exp16$year_implemented<=2016, 1, 0)
exp16$exp[is.na(exp16$year_implemented)] <- 0
exp16$year <- 2016
exp16 <- exp16[c("State","year","exp")]
#2018
exp18 <- exp
exp18$exp <- ifelse(exp18$year_implemented>=2014 & exp18$year_implemented<=2018, 1, 0)
exp18$exp[is.na(exp18$year_implemented)] <- 0
exp18$year <- 2018
exp18 <- exp18[c("State","year","exp")]
#2020
exp20 <- exp
exp20$exp <- ifelse(exp20$year_implemented>=2014 & exp20$year_implemented<=2020, 1, 0)
exp20$exp[is.na(exp20$year_implemented)] <- 0
exp20$year <- 2020
exp20 <- exp20[c("State","year","exp")]
#2012
exp12 <- exp
exp12$exp <- 0
exp12$year <- 2012
exp12 <- exp12[c("State","year","exp")]
#2010
exp10 <- exp
exp10$exp <- 0
exp10$year <- 2010
exp10 <- exp10[c("State","year","exp")]
#merging together
rm(exp)
exp <- rbind.data.frame(exp10,exp12,exp14,exp16,exp18,exp20)
kmerg.exp <- merge(kmerg, exp, by = c("State","year"))
# getting geo data for plots
data <- us_map("states")
names(kmerg.exp)[1] <- "full"
medplot <- merge(data,kmerg.exp, by = "full")
# convert var type for plots
medplot$exp <- as.character(medplot$exp)
# read in voting burden measures
vburd <- read.csv("burden_data/electburden_measures.csv")
# merge w geo data
names(vburd)[1] <- "abbr"
vburdplot <- merge(data,vburd, by = "abbr")

##### PLOTS #####
#### Figure 2 Medicaid burden over time ####
ggplot(data = medplot) + facet_wrap(~year, ncol=2) +
  geom_sf_pattern(aes(fill = negsum, pattern = exp), linewidth=.5, color="black") +
  theme_map() + theme(legend.position="top", legend.justification = "center") +
  scale_fill_continuous(name="Medicaid burden index", low="grey90", high="grey10", breaks = seq(0,1,.1), labels = dropLeadingZero) +
  scale_pattern_manual(values = c("0" = 'none',"1" = 'stripe'),name="Expanded Medicaid?", labels=c("No","Yes") ) +
  theme(legend.text = element_text(size=9),
        legend.title = element_text(size=10),
        strip.text.x = element_text(size=11, face="bold"),
        strip.background = element_rect(colour="white", fill="white"))
ggsave("fig2_medburden_allyears.png", width=7, height=9)

#### Figure 3 Average Med burden over time #### 
mean <- aggregate(kmerg$negsum, list(kmerg$year), FUN=mean)
names(mean) <- c("year","mean")

ggplot(mean, aes(x=year, y=mean)) + 
  geom_line() + geom_point() +
  theme_bw() + scale_x_continuous(breaks = seq(2010,2020,2)) +
  scale_y_continuous(breaks=seq(.2,.8,.1), limits = c(.2,.8)) +
  theme(panel.grid.major = element_blank(),
        legend.position = "none", 
        axis.text.x = element_text(colour="black",size="11"),
        axis.text.y = element_text(colour="black",size="11"),
        legend.text = element_text(size="12"),
        legend.title = element_text(size="12"),
        axis.title = element_text(size="12", face = "italic")) +
  ylab("Medicaid burden index") + xlab("Year") +
  geom_vline(xintercept = 2014, colour = "grey60", linetype = 2)

ggsave("fig3_index_over_time.png",width=4.5, height=4)

#### Figure 4 Voter registration burden over time ####
ggplot(data = vburdplot) + facet_wrap(~year, ncol=2) +
  geom_sf_pattern(aes(fill = regidx), linewidth=.5, color="black", pattern="none") +
  theme_map() + theme(legend.position="top", legend.justification = "center") +
  scale_fill_continuous(name="Registration burden index", low="grey90", high="grey10", breaks = seq(0,1,.2), labels = dropLeadingZero) +
  theme(legend.text = element_text(size=9),
        legend.title = element_text(size=10),
        strip.text.x = element_text(size=11, face="bold"),
        strip.background = element_rect(colour="white", fill="white"))
ggsave("fig4_regburden_allyears.png", width=7, height=9)

#### Figure 5 Voter turnout burden over time ####
ggplot(data = vburdplot) + facet_wrap(~year, ncol=2) +
  geom_sf_pattern(aes(fill = turnidx), linewidth=.5, color="black", pattern="none") +
  theme_map() + theme(legend.position="top", legend.justification = "center") +
  scale_fill_continuous(name="Turnout burden index", low="grey90", high="grey10", breaks = seq(0,1,.2), labels = dropLeadingZero) +
  theme(legend.text = element_text(size=9),
        legend.title = element_text(size=10),
        strip.text.x = element_text(size=11, face="bold"),
        strip.background = element_rect(colour="white", fill="white"))
ggsave("fig5_turnburden_allyears.png", width=7, height=9)

#cleaning up R environment
rm(list=ls()[!ls() %in% c("burd100")])

#### Figure 6 Coefficient plot ####
# load in results
load("mods_output.RData")
#cleaning up R environment
rm(list=ls()[!ls() %in% c("turn.bothidx.mod.100","turn.bothidx_full.mod.100","burd100")])
# load in results
load("bootstraps_mods_output.RData")
# keep indices only
t1m <- tidy(turn.bothidx.mod.100)[3,1:3]
t2m <- tidy(turn.bothidx_full.mod.100)[3,1:3]
t3m <- t1m
t3m[1,2] <- mod.tot.nocovar$coef$b
t3m[1,3] <- mod.tot.nocovar$results$Effects[2]
t4m <- t2m
t4m[1,2] <- mod.tot$coef$b
t4m[1,3] <- mod.tot$results$Effects[2]
# add model labels
t1m$model <- "Classic TWFE:\n burden + expansion only"
t2m$model <- "Classic TWFE:\n with covariates"
t3m$model <- "de Chaise & D'Hault 2024:\n burden + expansion only"
t4m$model <- "de Chaise & D'Hault 2024:\n with covariates"
# create plot data
t_comb <- rbind.data.frame(t1m,t2m,t3m,t4m)
plotdata <- t_comb
plotdata$term <- plotdata$model
# plot
dwplot(plotdata, dodge_size = .5,
       vline = geom_vline(xintercept = 0, colour = "grey60", linetype = 2)) +
  xlab("Coefficient estimate") + theme_bw() +
  theme(axis.title.x = element_text(colour="black",size = 9), axis.text.y = element_text(colour="black",size = 9), legend.position = "none") +
  scale_color_manual(values = rep("black",4)) 
ggsave("fig6_coefplot_indices.png", width=5, height=2.5)

#cleaning up R environment
rm(list=ls()[!ls() %in% c("burd100")])

#### Figure A1 pooled scatterplot #### 
ggplot(burd100,
       aes(y=turn,x=negsum)) + geom_jitter(alpha = 0.04) +geom_smooth(method="lm", color="black", linewidth=1.5) + 
  theme_bw() +
  theme(axis.text.x = element_text(colour="black",size="11"),
        axis.text.y = element_text(colour="black",size="11"),
        axis.title = element_text(size="12", face = "italic")) + ylab("County-level turnout") + xlab("Medicaid burden index")
ggsave("figA1_scatterplot.png",width=4.5, height=4)
